Simulation of stochastic quantum systems using polynomial chaos expansions 
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We present an approach to the simulation of quantum systems driven by classical stochastic processes 
that is based on the polynomial chaos expansion, a well-known technique in the field of uncertainty 
quantification. The polynomial chaos expansion represents the system density matrix as a series of 
orthogonal polynomials in the principle components of the stochastic process and yields a sparsely 
coupled hierarchy of linear differential equations. We provide practical heuristics for truncating 
this expansion based on results from time-dependent perturbation theory and demonstrate, via an 
experimentally relevant one-qubit numerical example, that our technique can be significantly more 
computationally efficient than Monte Carlo simulation. 



Introduction - Quantitative understanding of the dynam- 
ics of open quantum systems is critically important to 
many contemporary physics experiments pQ. While the 
equations of motion for open systems models are often 
simple to formulate, in only a few special cases may they 
be solved analytically, and numerical studies are often 
limited to small systems. However, in the absence of 
quantum back action, i.e., when the environmental dy- 
namics are independent of the system state, fully quan- 
tum open systems may be well approximated by semi- 
classical stochastic driving, whereupon the environment 
interaction operators are replaced by classical stochastic 
processes. For example, the coherence decay of diamond 
nitrogen-vacancy (NV) centers in the presence of dilute 
paramagnetic defects may be modeled very well by as- 
suming that paramagnetic defects in the lattice produce 
a classical, fluctuating Overhauser field which dephases 
the NV center [2 , 3 . Expensive numerical studies model- 
ing the full quantum environment are then only required 
only to determine the statistical properties of this effec- 
tive field. The resulting stochastic models are often suf- 
ficient to compute any desired system observables. How- 
ever, these reduced models exchange quantum degrees of 
freedom for stochastic ones that may also require large, 
but significantly reduced, computational overhead. Ex- 
pectation values of the system observables may then, in 
principle, be computed by averaging over the stochastic 
degrees of freedom in a manner that is consistent with 
the statistics of the stochastic process. In practice, how- 
ever, such an average is often difficult to calculate. Monte 
Carlo (MC) methods [1] approximate this average by gen- 
crating many sample noise trajectories, integrating the 
Schrodinger equation for each trajectory, and averaging 
the resulting density operators. However, MC can be no- 
toriously slow to converge, making it impractical for ap- 
plications requiring iterative numerical calculations, such 
as optimal control [5, 6J. Pcrturbative master equations, 
on the other hand, are often either computationally in- 
expensive and inaccurate, or expensive and accurate, de- 
pending on the approximations made in their derivation 
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In this work, we present an alternative approach to per- 
forming the stochastic average based on a class of tech- 
niques used widely in classical uncertainty quantification. 
Known as the polynomial chaos expansion (PCE) [7J[S], 
this method leverages properties of orthogonal polynomi- 
als to yield a converging sequence of approximate evolu- 
tion equations for a quantum system undergoing stochas- 
tic driving without resorting to MC methods. While we 
restrict our discussion to quantum systems driven by clas- 
sical Gaussian stochastic processes, we make no assump- 
tions of weak coupling nor do we restrict the form of the 
noise correlation function. Furthermore, we show that 
the linearity of the Schrodinger equation makes quantum 
systems particularly well suited to the PCE approach, as 
the stochastic dynamics may be expressed in terms of a 
sparsely-coupled system of differential equations. 

We begin this article with a derivation of Karhunen- 
Loeve decomposition, which expresses correlated, clas- 
sical stochastic processes as an easily truncated sum of 
deterministic functions weighted by uncorrelated random 
variables. We proceed to use this decomposition to de- 
rive the PCE as applied to stochastic quantum systems, 
yielding a sparsely-coupled system of Schrodinger-like 
equations. We conclude with a discussion and numerical 
simulation of the convergence properties of this method, 
benchmarking our results against Monte Carlo simula- 
tions. 

Model - We consider a quantum two-level system coupled 
linearly to a classical stochastic process, O, and described 
by the Hamiltonian, H(t;ti(t)) = H Q (t) +ti(t)V. Switch- 
ing to a rotating frame with respect to Hq, we obtain: 

H(t; n(t)) = n(t)u (tfvu = si(t)v(t), (1) 

where Uq = Texp(— i f Q Ho(s)ds) and T is the Dyson 

time-ordering operator. Hamiltonians of this form are 
quite common, and restriction to this minimal form sim- 
plifies the following derivations. Generalizations to mul- 
tiple or more complicated dependence on the stochastic 
process require only straightforward modifications to the 
following procedure. 



2 



We restrict our discussion here to stochastic processes, 
fl, which are mean-zero, Gaussian, and stationary [9]. By 
Wick's theorem [10], such processes may be completely 
described in terms of their two-point correlation func- 
tions, C(t\,t2) = (Sl(£i)Sl(£ 2 )) £1 . In this article, we use 
the notation (/(0(t))) n to signify t ne expectation value 
of the function / with respect to the process, f2. 

The state of the system when conditioned on a spe- 
cific realization of the stochastic process, p(t; {fl(t)}) , 
will evolve according to the Schrodinger-von Neumann 
equation: 



Eigenvalues 



Eigenfunctions 



f rfp(*; {"(*)» 

dt 



n(t)v(tr P (f,{m}), 



(2) 



where we have used the superoperator adjoint notation 
A x B = \A,B\. At a time r, the state of the system, av- 
eraged over the stochastic process, is given by the formal 
expression 

p(r) = (U(t; m)})p(0)U(r; {0(i)})) n , (3) 

where J7(r;{0(t)}) denotes the unitary operator gener- 
ated by Eq. ^ with a specific realization of fl. The ob- 
jective of this work is to demonstrate that this stochastic 
average may be performed in a computationally efficient 
manner using the PCE. 

Karhunen-Loeve Expansion - If the stochastic process, 
57, is white, the average in Eq. (§ may be taken locally 
in time, and the system dynamics may be described ex- 
actly by a Lindblad master equation [T|. However, the 
presence of non- vanishing time correlations greatly com- 
plicates the calculation of the stochastic average. To sim- 
plify this calculation, we employ the Karhunen-Loeve ex- 
pansion (KLE) ,8J, which expresses a continuous, corre- 
lated process in terms of a discrete sum of deterministic 
functions weighted by uncorrelated random-variables: 



(4) 



Here, £ n € A/"(0, 1) are independent and identically 
distributed (iid) random variables drawn from a unit- 
variance, zero-mean Gaussian distribution, while A„ 
and g n (t) are, respectively, the eigenvalues and L 2 - 
orthonormal eigenfunctions of the Fredholm equation [8] : 



C(ti,t2)g n (t2)dt 2 = \ n g n (ti)- 



(5) 



Here, the correlation function acts as a symmetric, posi- 
tive scmi-dcfinitc integral kernel, so Mercer's theorem [8 
implies that the eigenvalues, A„, are discrete and non- 
negative. Non-negativity is ensured because the correla- 
tion functions of stationary processes are positive semi- 
definite (by Bochner's theorem [H]), while discreteness 
is guaranteed by the finite upper limit on the integral, 
Eq. f5|. 




Figure 1. (color online) The eigenvalues, X n (left, gray x's); 
cumulative transition rates, r„, (left, solid circles); and eigen- 
functions, g n (t) (right), of the Fredholm equation, Eq. Q5B, 
for Orenstein-Uhlenbeck noise, where C(t) = exp(— \t\ /tc), 
calculated with final time r = 1 and Hamiltonian H(t) = 
Bo x + Q(t)a z . The top row corresponds to noise correlation 
time t c = 0.1 and magnetic field B — 5, while the bottom row 
for t c = 10 and magnetic field B — 20. The stochastic modes 
with the three largest transition rates are color coordinated 
with the eigenfunction plots (in descending transition rates: 
red dot and long dashed line; green dot and short dashed line; 
blue dot and solid blue line). For long correlation times, com- 
pared to r, only one mode is dominant, while more modes are 
increasingly important for shorter correlation times. 



Interestingly, when the final time is much greater 
than the correlation time of the stochastic process, i.e., 
t 3> t c , the Wiener-Khinchin theorem |10j implies that 
the eigenvalue spectrum becomes continuous and equal 
to the noise power spectral density, i.e., A w = S(oj), 
while the eigenfunctions take the form g u {t) cx cos(cjf). 
Taken to the extreme white-noise limit, where t c — > 0, 
all eigenvalues are equal. In the regime where r — > r c , 
the correlation function is approximately constant over 
the integration window, and the Fredholm equation pos- 
sesses only a single nonzero eigenvalue, with eigenfunc- 
tion g\(t) oc 1. In this case, the stochastic process may 
be well approximated as a random variable which is con- 
stant over t g [0, t]. Figure [l] illustrates these two limits 
by plotting the eigensystem of the Fredholm equation for 
Orenstein-Uhlenbeck [10] noise with two different decay 
parameters. 

When the correlation time is relatively long compared 
to the evolution time, but not infinite, e.g., when r < r c , 
the expansion Eq. Q is dominated by only a few terms 
corresponding to the largest S eigenvalues, and so may 
be truncated with minimal error. However, truncating 
the KLE based only on the eigenvalues ignores the im- 
pact that higher frequency modes could have on the dy- 
namics, e.g.,, resonance. We therefore propose a more 
physically motivated truncation criterion based on results 
from time-dependent perturbation theory |12j . For static 
Hamiltonian terms Hq and V in the Schrodinger picture, 
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the rate at which any given mode, g n (t), could cause a 
transition between the jth and fcth eigenstates of Hq is 
given by 



Y° = 



1 



U\v\k) i , 

o 



- Ek)t y/K t g n {t)dt 



(6) 



where i?ob) = -Ejb')- Summing over these eigenstates 
provides a measure of the degree to which a given mode 
will impact the evolution of the system, i.e., the cumu- 
lative transition rate, T n = ^ . fc . In addition to the 
eigenvalues, A„, and eigenfunctions, g n (t), of the Fred- 
holm equation, Eq. ([5]), cumulative transition rates are 
also included in Figure [T] Thus, we approximate the ex- 
pansion Eq. Q by keeping only those modes correspond- 
ing to the S largest transition rates, T„; we shall refer to 
S as the stochastic dimension. With this approximation, 
Eq. |2| becomes 



dt 



s 

E 

n=l 



X n g n mnV(t) x p(t;0- 



(7) 



We emphasize that this truncation strategy differs from 
that usually taken in standard uncertainty quantification 
literature [5] , wherein the truncation is based on the mag- 
nitude of the eigenvalues alone, and without considera- 
tion to the potential impact of a given stochastic mode 
on the system dynamics. 

Expansion in orthogonal polynomials - At the final time, 
t, the state of the system may be considered as a com- 
plicated function of the S uncorrelated random variables 
from Eq. Q, i.e., p(r) = p(t;|), as expressed in Eq. 
As such, we may expand this function in a complete basis 
of orthogonal polynomials, $„(£): 



p(*;£) = 5>n(*)$»(l') I 



(8) 



n=0 



yielding the (untruncated) PCE. Here, 4> n (t) are the 
time-dependent expansion coefficients that represent our 
new dynamical variables. The polynomials should be or- 
thogonal under the measure, = exp(— £ 2 /2)/V27rd£, 
which is derived from the stationary distribution of the 
random variables, which are drawn from A/"(0, 1). Multi- 
variate Hermite polynomials are the natural choice: 
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where we now use the multi- index vector n 6 Z> . This 
expansion, Eq. Q, may be truncated, keeping only terms 
for which Hnl^ = J2j n j — P-> where P is the PCE or- 
der. Inserting the truncated expansion into the evolution 
equation, Eq. Q, we have 



||k|| 1= 



dt 



s 

E 



P 

E< 

l|l|li=0 



Exploiting the orthogonality of the Hermite polynomials, 
we may compute the evolution equation for the coeffi- 
cients of the expansion of Eq. (l8|): 



.d4> m {t) 
dt 



s p 

= E E \^9n(t)V(t) x Mt)G mnh (9) 

n=l||l|| 1=0 



where we have defined the Galerkin projection: 



*m(0£„*l(0 



(10) 



These projection terms may be computed explicitly using 
the Hermite polynomial orthogonality relations: 

H.e n (x)He m (x)e~ x ^ 2 dx — y/2xnl6 min , 

> 

and the three term recurrence relation: 

£He„(£) = He„ +1 (£) + uRc^d). 

Taken together, these lead to a much-simplified expres- 
sion for the Galerkin projection: 



G mn \ — ((m„ + l)£m„+i,i n 



-1,0 H ^m 3 l, 



Owing to the presence of the Kronecker delta functions 
in the Galerkin projection, the PCE results in a sparsely- 
coupled hierarchy of linear differential equations. The 
choice of both the stochastic dimension S and the PCE 
order P determine the number of equations N in the 
hierarchy through a simple combinatorial argument [8]: 



N 



p 

E 

m=0 



S + m — 1 
5-1 



(S + P)\ 

s\p\ ■ 



(11) 



This scaling, known colloquially as the curse of dimen- 
sionality, limits practical applications to those situations 
in which i) the noise correlation time is long, resulting 
in low stochastic dimension and ii) the noise is weak, 
so that the PCE converges quickly. More sophisticated 
truncation procedures may reduce the hierarchy depth, 
however, this remains an area of active research. 
Convergence of the PCE - The convergence properties 
of the coupled evolution equations, Eq. depend crit- 
ically on the distribution of the cumulative transition 
rates, r„. Specifically, noise modes associated with large 
transition rates will couple strongly to the system and the 
PCE must be truncated at high order in those variables 
to faithfully represent the system dynamics. For exam- 
ple, modes for which T n r > 1 will, on average, induce 
at least one transition over the course of the evolution. 
Accurately capturing these dynamics would require such 
modes to be considered at high PCE order. 
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Figure 2. (Color online) A comparison of the predicted time- 
dependent coherence, (a x (t)), for the Hamiltonian given in 
Eq. (12 I, with correlation function C(t) — 9 exp(— \t\ /10). 



Solid, colored lines correspond to a stochastic dimension S = 
1, while dashed lines correspond to S = 3. Grey line is the 
Monte Carlo result drawn with width equal to the standard 
error in the estimate. Because of the long correlation time, 
only one eigenvalue of the KLE is dominant, so it is more 
efficient to keep the stochastic dimension small and increase 
the PCE order. 



We consider explicitly the stochastic dynamics of a 
driven quantum two-level system, or qubit, coupled to 
a classically fluctuating dephasing process: 



H(t) =a x + n(t)a z 



(12) 



where a x and a z are Pauli matrices. Such a model de- 
scribes, for example, Rabi oscillations in the presence of 
dephasing noise [H [13] , and is particularly relevant for 
NV centers in diamond [2], for example. In the absence 
of the drift term a x , the dephasing dynamics are exactly 
solvable for any stationary Gaussian process, Q. How- 
ever, when this term is included, the Hamiltonian does 
not commute with itself at different times, and the sys- 
tem is no longer analytically integrable. To illustrate 
our method, we choose the stochastic process, fi, to be 
Gaussian Orenstein-Uhlenbeck type [T3] , with correlation 
function of the form C(t) = a 2 exp(— \t\ /r c ); the cou- 
pling parameter, a, and the correlation time, r c , will be 
specified later. We specify the initial state of the system 
as |o-+), where a x \af) = ±\af), and we compute the 
time-dependent coherence, (a x (t)). By tuning the noise 
correlation time and the coupling parameter, this model 
can explore the convergence of our method with respect 
to PCE order P and stochastic dimension S. 

To benchmark the performance of our PCE approach, 
we compare against MC simulations. MC algorithms ap- 
proach the problem of computing the stochastic average 
in Eq. ([3| by generating a sufficiently large number of 
statistically consistent realizations of the stochastic pro- 
cess, Q, evolving the system with each of the realizations, 
and averaging the final-time density matrices. 

As shown in Fig. [2j our PCE method is capable of 



reproducing the results of MC simulations with high ac- 
curacy, significantly faster than MC. For the example 
chosen, Monte Carlo required approximately 4000 itera- 
tions for convergence, while the most accurate PCE re- 
sults report here (P = 9 and S = 3) required the solution 
of only 220 coupled equations and ran approximately 20 
times faster than the MC simulation. Note that because 
r c /r = 10 in our simulation, only one eigenvalue of the 
KLE is dominant. In this regime, it is more efficient to 
keep the stochastic dimension small (S < 3) and increase 
the PCE order for improved accuracy. As the order in- 
creases from P = 1 to P = 9, the accuracy of the PCE 
coherence increases as a function of time, compared to 
the converged MC result. 

Discussion - Our PCE method demonstrates the ability 
to rapidly and accurately propagate stochastic quantum 
systems. It outperforms MC simulations in computa- 
tional efficiency, and has the potential to become an im- 
portant tool in the study of noisy quantum systems. An 
area in which we expect the PCE to be particularly useful 
is in the realm of optimal control (OC). The high compu- 
tational cost of MC simulations severely limits its use in 
sequential optimal control simulations. However, PCEs 
can be both accurate and fast, and they may be easily in- 
corporated as part of a surrogate dynamical model in OC 
simulations. In this work, the PCE has been formulated 
to propagate a particular state, however, the equation of 
motion for the complete dynamical map takes a similar 
form and may also be adapted easily to PCE methods. 
Implementation of state-to-state and dynamical map OC 
will appear in forthcoming work. 

An interesting comparison can be made between the 
PCEs as presented here and another expansion based 
on orthogonal polynomials: Kubo's hierarchy equations 
of motion (HEOM) [15], and their generalization to all 
diffusive processes, the DHEOM [T5]. Application of 
the DHEOM /HEOM demands that the noise be diffu- 
sive and have exponentially decaying correlation func- 
tion, and proceeds by diagonalizing the noise generating 
functional using of orthogonal polynomials [16] . Inter- 
estingly, though the HEOM and PCE approaches each 
yields a sparsely-coupled hierarchy of differential equa- 
tions based on expansions in orthogonal polynomials, 
they perform well in exactly opposite limits: the HEOM 
method converges quickly for noise with a short corre- 
lation time, while the stochastic dimension of our PCE 
method converges quickly for noise with a long correla- 
tion time. For systems coupled to multiple, uncorrelated 
noise sources, it may be computationally advantageous 
to apply different methods for each source: HEOM for 
noise with short correlation times, PCEs for noise with 
long correlation times. 

Additionally, we have assumed a linear coupling be- 
tween the system and the stochastic process in Eq. ([!]). 
While such a coupling is common |17j , nonlinear interac- 
tions are possible, which may increase the coupling den- 
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sity of the differential equations in Eq. Q by modifying 
the form of the Galerkin projection of Eq. ( 10 ). Strongly 



nonlinear interactions will yield densely coupled systems 
of equations, which will increase the computational cost 
of this method. 

Despite the obvious utility of our PCE approach for 
simulating stochastic quantum systems, it does have lim- 
itations. Principal among these is the uncontrolled ap- 
proximation of the stochastic average, so that the error 
must be estimated by increasing the PCE order and/or 
the stochastic dimension of the expansion until conver- 
gence is seen. Furthermore, as indicated in Eq. (Ill, 



the number of equations to be solved grows combinato- 
rially with both the stochastic dimension and the PCE 
order; when these are large, the method of Galerkin pro- 
jections becomes computationally infeasibile. Intermedi- 
ate between MC and the PCE approach presented here 
is a non-intrusive formulation of the PCE, so called be- 
cause its implementation requires only a forward solver 
for the equations of motion, while the intrusive method 
presented here requires an explicit solver to propagate 
the coupled equations resulting from the Galerkin projec- 
tions. Non-intrusive spectral methods approximate the 
stochastic average of Eq. (|3| by performing a KLE, and 
using sparse quadrature techniques to perform the aver- 
age. We plan to implement such techniques in the near 
future. 
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